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DECUS Program Library Write-up 


DECUS NO. 8-483 


GRFIT, A Simple Least Squares Routine 
by R. C. Gross 

In the environment in which many small computers are found, 
there is often apparent a need for a simple statistical routine 
such as found in GRFIT but not the facilities to write or debug 
the program. While GRFIT is essentially similar to many other 
available least squares type programs, this is the first time, to 
my knowledge, that a similar program has been made available for 
PDP-8 machines in both Fortran and Focal.* In addition, the Fortran 
program has been written both as a main program and as a subroutine. 

In FORTRAN the program requires a PDP-8 with 8K of memory and 
one large storage device, either the single DECTAPE or a DF-32 disk 
so that the facility supports the PS-8 operating system. 

In FOCAL the program requires only the minimum configuration 
PDP-8 with 4K of memory. 

The FORTRAN program as dimensioned allows an additional 16 
pages of memory free for other material. It can, of course, be 
considerably reduced in its dimension statement. For users really 
short of space, the loop following the blank comment line ending 
at statement 12 can be removed; this removes the actual calculation 
of both the calculated values of Y and the differences between the 
calculated and the entered values of Y as well as cutting the 
dimension requirements in half. The calculated slopes and inter¬ 
cepts as well as most of the error information are still available 

* See FOCAL8-209 
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to the user. The forceoff statement immediately preceding statement 
30 in the main program can be changed to any value not necessary as 
a potential X value if 999 is not satisfactory. 

The FOCAL program accepts a much more limited array size, but 
it will easily accept the 8 to 12 pairs of points that are generally 
used with such a program. If 8-K FOCAL were available the array size 
would, of course, be considerably expanded. 

I emphasize that the programs are simple, the sort that are 
available at every sophisticated installation, but might be of use 
in the somewhat unusual environment in which many PDP-8 computers 
are installed. 

Variables; 

X - Array of input x values 
Y - Array of input y values 
YC - Array of calculated y values 

DIFF - Array of calculated differences between Y and YC 

NV - Number of points 

B - Calculated slope 

A - Calculated intercept 

R - Correlation coefficient 

SUM - Sum of the squares of the residuals 

EA - Standard error in intercept 

EB - Standard error in slope 

Note that Arrays of X and Y, as well as a value for NV, must 
be passed to the subroutine, while the remainder of the variables 
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are calculated and returned to the main program by the subroutine. 
The dimensioning of the main program must account for the returned 
arrays YC and DIFF. 

Cautionary Notes: 

In any statistical program the values of the calculated errors 
must be examined carefully. In the case of the limitations of 12 
bit words the accuracy of the results must always be examined care¬ 
fully, especially when small differences are involved. 

Values of YC and DIFF can be produced in the main program with¬ 
out the accompanying dimensioning of them as arrays if they are 
printed immediately after their calculation and never stored. 
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DIMENSION X C1 35 ), Y C 1 35 ) , YCC 1 35 ) , 1)1 KKC135 ) 

I N= 1 
I OUT:l 

WRITE(I OUT,105) 

WRITE(I OUT,107) 

SX = 0. 

SY = 0. 

SY2 = 0. 

N V = 0 
SXY = 0. 

SX2 = 0. 

DO 30 J:1,200 

NV = NV + 1 

READCIN,100)X(J),Y(J) 

IF CXCJ)-999.)30,40,30 
30 CONTINUE 
40 NV = NV-1 

TV:FLOATC NV) 

DO 11 I:1,NV 
SX : SX + X(I) 

SY : SY + YCI) 

SXY : SXY + X(I) * Y(I) 

SY2 : SY2 + YCI) * YCI) 

11 SX2 : SX2 + X(I) * X(I) 

SSX : SX * SX 

B:(TV*SXY-SX*SY)/(TV*SX2-SSX) 

A : C SY-B*SX)/TV 

R:CTV*SXY-SX*SY)/CSQRTCTV*SX2-SX*SX)*SQRTCTV*SY2-SY*SY)) 
RR:SQRT(CSY2-A*SY-B*SXY)/CTV-2.) ) 
EA:RR*SQRTCSX2/CTV*SX2-SX*SX)) 

£B:RR*SQRT(TV/(TV*SX2-SX*SX)) 

SUM : 0. 

C 

DO 12 1:1,NV 
YCCI) : A+ B*XCI) 

DIFFCI) : YCCI) - YCI) 

SUM = SUM + DIFFCI) * DIFFCI) 

12 CONTINUE 
WRITECIOUT,106) 

DO 50 1:1,NV 

WRITEClOUT,210)XCI),YCI),YCCI),DIFFCI) 

50 CONTINUE 

105 FORMATC44H INPUT FORMAT IS 2F10.4, X:999. TERMINATES ) 

WRITECIOUT,215)B,A 

WRITECIOUT,216)R,SUM 
WRITECIOUT,218)EB,EA 

218 FORMATC’ SLOPE ST ERROR: ’ ,E10.3, ’ INTERCEPT ST ERROR:’,E10.3) 

215 F0RMATC7H SLOPE:FI 0.4,11H INTERCEPT:F10.4) 

216 FORMATC13H CORR.COEFF :F10.4,18H SUM SQ.RESID.ERR:F10.4) 

106 FORMATC 43H XCI) YCI) YCCI) DIFFCI) ) 

100 FORMATC2F10.4) 

210 FORMATC4F10.4) 

107 FORMATC’ X Y’) 

CALL EXIT 

END 
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Fortran Subroutine Listing 


SUBROUTINE GRFI TCX*Y>NWYC.» DIFF* B» A# K# SUM> EA# EB) 

DI MEN SION XC1)*YC1)*YCC1)# BIFFC 1 ) 

SX = 0. 

SY = 0. 

SY2 = 0. 

NV = 0 
SXY = 0. 

SX2 = 0. 

IV = FLOAT CNV ) 

DO 11 I = 1« N V 
SX = SX + ACI) 

SY = SY + Y(I) 

SXY = SXY + A(I) * YCI) 

SY2 = SY2 ♦ YCI) * YCI) 

SA2 = SX2 + XCI) * AC I) 

SSX = SX * SX 

B=CTV*SXY-SX*SY)/C1V+SX2-SSA) 

A = C SY-B+SX)/1V 

K= C I V+SXY-SX+S Y ) / C SQR1 C 1 V*iA2-SX*SX)*S(df<TC TV*SY2-SY*SY)) 

SUM = 0. 

HH=SQR1 C CSY2-A*SY-B*SX,Y)/,CT,V.-2- ) ). 

EA=Rk* SQkT C SX2/C 1 V/*SX2-SA*SA) ) 
EB=KK*SQHTC1V/CTV*SX2-SX*SX)) 

DO 12 I = 1#NV 
YCCI) = A+ B*XCI) 

DIFFCI) = YCCI) - YCI) 

SUM = SUM + DIFFCI) * DIFFCI) 

CON TINUE 
RETURN 
CALL EXIT 
END 
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Typical Data Set - Fortran 


Input 

x 


format 


0. 

1« 


999. 


0 

1 

2 < 

3. 

4. 


XCI ) 
■ 0000 
> 0000 
0000 
0000 
0000 


IS 

Y 


2F10.4* X=999• 


TERMINATES 


0« 


1. 

2. 08 
2. 95 

4. 


Y ( I ) 
0.0000 
0000 
0800 
9500 
0000 


YC(I) 
0.0160 
01 10 
0060 
0010 
9960 


SLOPE= 0.9950 INTERCEPT 5 
CORR.COEFF = 0.9996 SUM 

SLOPE ST ERROR 5 0.168E-01 


DIFFCI) 
0. 0160 
01 10 
0740 
0510 
0040 


0 

-0 
0 

-0 
0. 0160 

SQ.RESID.ERR 5 


0 


INTERCEPT ST ERROR= 


. 0085 
0. 41 IE 
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